Открытый курс по машинному обучению. Сессия № 2

Автор материала: Data Science интерн Ciklum, студент магистерской программы CSDS UCU Виталий Радченко, программист-исследователь Mail.ru Group, старший преподаватель Факультета Компьютерных Наук ВШЭ Юрий Кашницкий. Материал распространяется на условиях лицензии Creative Commons CC BY-NC-SA 4.0. Можно использовать в любых целях (редактировать, поправлять и брать за основу), кроме коммерческих, но с обязательным упоминанием автора материала.

Домашнее задание №5

Случайный лес и логистическая регрессия в задачах кредитного скоринга и классификации отзывов к фильмам

Нашей главной задачей будет построение и настройка моделей для задач кредитного скоринга и анализа отзывов к фильмам. Заполните код в клетках (где написано "Ваш код здесь") и ответьте на вопросы в веб-форме.

Но для разминки решите первое задание.

Задание 1. В зале суда есть 7 присяжных, каждый из них по отдельности с вероятностью 80% может правильно определить, виновен подсудимый или нет. С какой вероятностью присяжные все вместе вынесут правильный вердикт, если решение принимается большинством голосов?

Варианты ответа:

  • 20.97%
  • 80.00%
  • 83.70%
  • 96.66%

Решение: поскольку большинство голосов – 4, тогда наше m = 4, N = 7, p = 0.8. Подставляем в формулу из статьи $$ \large \mu = \sum_{i=4}^{7}C_7^i0.8^i(1-0.8)^{7-i} $$ После подставления и проделывания всех операций получим ответ 96.66%


In [1]:
import math

def nCr(n,r):
    f = math.factorial
    return f(n) / f(r) / f(n - r)

p, N, m, s = 0.8, 7, 4, 0
for i in range(m, N+1):
    s += nCr(N, i) * p**i * (1 - p) ** (N - i)
print(s)


0.9666560000000001

Теперь перейдем непосредственно к машинному обучению.

Данные по кредитному скорингу представлены следующим образом:

Прогнозируемая переменная
  • SeriousDlqin2yrs – Человек имел долгие просрочки выплат платежей за 2 года; бинарный признак
Независимые признаки
  • age – Возраст заёмщика кредитных средств (число полных лет); тип – integer
  • NumberOfTime30-59DaysPastDueNotWorse – Количество раз, когда человек имел просрочку выплаты других кредитов более 30-59 дней (но не больше) в течение последних двух лет; тип - integer
  • DebtRatio – Ежемесячный отчисления на задолжености(кредиты,алименты и т.д.) / совокупный месячный доход percentage; тип – float
  • MonthlyIncome – Месячный доход в долларах; тип – float
  • NumberOfTimes90DaysLate – Количество раз, когда человек имел просрочку выплаты других кредитов более 90 дней; тип – integer
  • NumberOfTime60-89DaysPastDueNotWorse – Количество раз, когда человек имел просрочку выплаты других кредитов более 60-89 дней (но не больше) в течение последних двух лет; тип – integer
  • NumberOfDependents – Число человек в семье кредитозаёмщика; тип – integer

In [2]:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline

Напишем функцию, которая будет заменять значения NaN на медиану в каждом столбце таблицы.


In [3]:
def impute_nan_with_median(table):
    for col in table.columns:
        table[col]= table[col].fillna(table[col].median())
    return table

Считываем данные


In [4]:
data = pd.read_csv('../../data/credit_scoring_sample.csv', sep=";")
data.head()


Out[4]:
SeriousDlqin2yrs age NumberOfTime30-59DaysPastDueNotWorse DebtRatio NumberOfTimes90DaysLate NumberOfTime60-89DaysPastDueNotWorse MonthlyIncome NumberOfDependents
0 0 64 0 0.249908 0 0 8158.0 0.0
1 0 58 0 3870.000000 0 0 NaN 0.0
2 0 41 0 0.456127 0 0 6666.0 0.0
3 0 43 0 0.000190 0 0 10500.0 2.0
4 1 49 0 0.271820 0 0 400.0 0.0

Рассмотрим типы считанных данных


In [5]:
data.dtypes


Out[5]:
SeriousDlqin2yrs                          int64
age                                       int64
NumberOfTime30-59DaysPastDueNotWorse      int64
DebtRatio                               float64
NumberOfTimes90DaysLate                   int64
NumberOfTime60-89DaysPastDueNotWorse      int64
MonthlyIncome                           float64
NumberOfDependents                      float64
dtype: object

Посмотрим на распределение классов в зависимой переменной


In [6]:
ax = data['SeriousDlqin2yrs'].hist(orientation='horizontal', color='red')
ax.set_xlabel("number_of_observations")
ax.set_ylabel("unique_value")
ax.set_title("Target distribution")

print('Distribution of target:')
data['SeriousDlqin2yrs'].value_counts() / data.shape[0]


Distribution of target:
Out[6]:
0    0.777511
1    0.222489
Name: SeriousDlqin2yrs, dtype: float64

Выберем названия всех признаков, кроме прогнозируемого


In [7]:
independent_columns_names = data.columns.values
independent_columns_names = [x for x in data if x != 'SeriousDlqin2yrs']
independent_columns_names


Out[7]:
['age',
 'NumberOfTime30-59DaysPastDueNotWorse',
 'DebtRatio',
 'NumberOfTimes90DaysLate',
 'NumberOfTime60-89DaysPastDueNotWorse',
 'MonthlyIncome',
 'NumberOfDependents']

Применяем функцию, заменяющую все значения NaN на медианное значение соответствующего столбца.


In [8]:
table = impute_nan_with_median(data)

Разделяем целевой признак и все остальные – получаем обучающую выборку.


In [9]:
X = table[independent_columns_names]
y = table['SeriousDlqin2yrs']

Бутстрэп

Задание 2. Сделайте интервальную оценку (на основе бутстрэпа) среднего дохода (MonthlyIncome) клиентов, просрочивших выплату кредита, и отдельно – для вовремя заплативших. Стройте 90% доверительный интервал. Найдите разницу между нижней границей полученного интервала для не просрочивших кредит и верхней границей – для просрочивших. То есть вас просят построить 90%-ые интервалы для дохода "хороших" клиентов $[good\_income\_lower, good\_income\_upper]$ и для "плохих" – $[bad\_income\_lower, bad\_income\_upper]$ и найти разницу $good\_income\_lower - bad\_income\_upper$.

Используйте пример из статьи. Поставьте np.random.seed(17). Округлите ответ до целых.

Варианты ответа:

  • 344
  • 424
  • 584
  • 654

In [10]:
def get_bootstrap_samples(data, n_samples, seed=0):
    # функция для генерации подвыборок с помощью бутстрэпа
    np.random.seed(seed)
    indices = np.random.randint(0, len(data), (n_samples, len(data)))
    samples = data[indices]
    return samples

def stat_intervals(stat, alpha):
    # функция для интервальной оценки
    boundaries = np.percentile(stat, [100 * alpha / 2., 100 * (1 - alpha / 2.)])
    return boundaries

# сохранение в отдельные numpy массивы данных по просрочке
churn = data[data['SeriousDlqin2yrs'] == 1]['MonthlyIncome'].values
not_churn = data[data['SeriousDlqin2yrs'] == 0]['MonthlyIncome'].values

# генерируем выборки с помощью бутстрэра и сразу считаем по каждой из них среднее
churn_mean_scores = [np.mean(sample) 
                     for sample in get_bootstrap_samples(churn, 1000, seed=17)]
not_churn_mean_scores = [np.mean(sample) 
                         for sample in get_bootstrap_samples(not_churn, 1000, seed=17)]

#  выводим интервальную оценку среднего
print("Mean interval",  stat_intervals(churn_mean_scores, 0.1))
print("Mean interval",  stat_intervals(not_churn_mean_scores, 0.1))
print("Difference is", stat_intervals(not_churn_mean_scores, 0.1)[0] - 
      stat_intervals(churn_mean_scores, 0.1)[1])


Mean interval [ 5462.17301516  5641.20421404]
Mean interval [ 6295.93237577  6505.35467934]
Difference is 654.728161731

Дерево решений, подбор гиперпараметров

Одной из основных метрик качества модели является площадь под ROC-кривой. Значения ROC-AUC лежат от 0 до 1. Чем ближе значение ROC-AUC к 1, тем качественнее происходит классификация моделью.

Найдите с помощью GridSearchCV гиперпараметры DecisionTreeClassifier, максимизирующие площадь под ROC-кривой.


In [11]:
from sklearn.tree import DecisionTreeClassifier
from sklearn.model_selection import GridSearchCV, StratifiedKFold

Используем модуль DecisionTreeClassifier для построения дерева решений. Из-за несбалансированности классов в целевом признаке добавляем параметр балансировки. Используем также параметр random_state=17 для воспроизводимости результатов.


In [12]:
dt = DecisionTreeClassifier(random_state=17, class_weight='balanced')

Перебирать будем вот такие значения гиперпараметров:


In [13]:
max_depth_values = [5, 6, 7, 8, 9]
max_features_values = [4, 5, 6, 7]
tree_params = {'max_depth': max_depth_values,
               'max_features': max_features_values}

Зафиксируем кросс-валидацию: стратифицированная, 5 разбиений с перемешиванием, не забываем про random_state.


In [14]:
skf = StratifiedKFold(n_splits=5, shuffle=True, random_state=17)

Задание 3. Сделайте GridSearch с метрикой ROC AUC по гиперпараметрам из словаря tree_params. Какое максимальное значение ROC AUC получилось (округлите до 2 знаков после разделителя)? Назовем кросс-валидацию устойчивой, если стандартное отклонение метрики качества на кросс-валидации меньше 1%. Получилась ли кросс-валидация устойчивой при оптимальных сочетаниях гиперпараметров (т.е. обеспечивающих максимум среднего значения ROC AUC на кросс-валидации)?

Варианты ответа:

  • 0.82, нет
  • 0.84, нет
  • 0.82, да
  • 0.84, да

In [15]:
dt_grid_search = GridSearchCV(dt, tree_params, n_jobs=-1, scoring ='roc_auc', cv=skf)
dt_grid_search.fit(X, y)


Out[15]:
GridSearchCV(cv=StratifiedKFold(n_splits=5, random_state=17, shuffle=True),
       error_score='raise',
       estimator=DecisionTreeClassifier(class_weight='balanced', criterion='gini',
            max_depth=None, max_features=None, max_leaf_nodes=None,
            min_impurity_decrease=0.0, min_impurity_split=None,
            min_samples_leaf=1, min_samples_split=2,
            min_weight_fraction_leaf=0.0, presort=False, random_state=17,
            splitter='best'),
       fit_params=None, iid=True, n_jobs=-1,
       param_grid={'max_depth': [5, 6, 7, 8, 9], 'max_features': [4, 5, 6, 7]},
       pre_dispatch='2*n_jobs', refit=True, return_train_score=True,
       scoring='roc_auc', verbose=0)

In [16]:
round(float(dt_grid_search.best_score_), 3)


Out[16]:
0.82

In [17]:
dt_grid_search.best_params_


Out[17]:
{'max_depth': 7, 'max_features': 6}

In [18]:
dt_grid_search.cv_results_["std_test_score"][np.argmax(dt_grid_search.cv_results_["mean_test_score"])]


Out[18]:
0.0027584835535929349

Простая реализация случайного леса

Задание 4. Реализуйте свой собственный случайный лес с помощью DecisionTreeClassifier с лучшими параметрами из прошлого задания. В нашем лесу будет 10 деревьев, предсказанные вероятности которых вам нужно усреднить.

Краткая спецификация:

  • Используйте основу ниже
  • В методе fit в цикле (i от 0 до n_estimators-1) фиксируйте seed, равный (random_state + i). Почему именно так – неважно, главное чтоб на каждой итерации seed был новый, при этом все значения можно было бы воспроизвести
  • Зафиксировав seed, выберите без замещения max_features признаков, сохраните список выбранных id признаков в self.feat_ids_by_tree
  • Также сделайте bootstrap-выборку (т.е. с замещением) из множества id объектов
  • Обучите дерево с теми же max_depth, max_features и random_state, что и у RandomForestClassifierCustom на выборке с нужным подмножеством объектов и признаков
  • Метод fit возвращает текущий экземпляр класса RandomForestClassifierCustom, то есть self
  • В методе predict_proba опять нужен цикл по всем деревьям. У тестовой выборки нужно взять те признаки, на которых соответсвующее дерево обучалось, и сделать прогноз вероятностей (predict_proba уже для дерева). Метод должен вернуть усреднение прогнозов по всем деревьям.

Проведите кросс-валидацию. Какое получилось среднее значение ROC AUC на кросс-валидации? Выберите самое близкое значение.

Варианты ответа:

  • 0.823
  • 0.833
  • 0.843
  • 0.853

In [19]:
from sklearn.base import BaseEstimator
from sklearn.model_selection import cross_val_score

class RandomForestClassifierCustom(BaseEstimator):
    def __init__(self, n_estimators=10, max_depth=10, max_features=10, 
                 random_state=17):
        self.n_estimators = n_estimators
        self.max_depth = max_depth
        self.max_features = max_features
        self.random_state = random_state
        
        self.trees = []
        self.feat_ids_by_tree = []
        
    def fit(self, X, y):
        for i in range(self.n_estimators):
            
            np.random.seed(i + self.random_state)
            
            feat_to_use_ids = np.random.choice(range(X.shape[1]), self.max_features, 
                                              replace=False)
            examples_to_use = list(set(np.random.choice(range(X.shape[0]), X.shape[0],
                                              replace=True)))
            
            self.feat_ids_by_tree.append(feat_to_use_ids)
            
            dt = DecisionTreeClassifier(class_weight='balanced',
                                        max_depth=self.max_depth, 
                                        max_features=self.max_features, 
                                        random_state = self.random_state)

            dt.fit(X[examples_to_use, :][:, feat_to_use_ids], y[examples_to_use])
            self.trees.append(dt)
        return self
    
    def predict_proba(self, X):
        predictions = []
        for i in range(self.n_estimators):
            feat_to_use_ids = self.feat_ids_by_tree[i]
            predictions.append(self.trees[i].predict_proba(X[:,feat_to_use_ids]))
        return np.mean(predictions, axis=0)

In [20]:
rf = RandomForestClassifierCustom(max_depth=7, max_features=6).fit(X.values, y.values)

In [21]:
cv_aucs = cross_val_score(RandomForestClassifierCustom(max_depth=7, max_features=6), 
                          X.values, y.values, scoring="roc_auc", cv=skf)
print("Средняя ROC AUC для собственного случайного леса:", np.mean(cv_aucs))


Средняя ROC AUC для собственного случайного леса: 0.831432621991

Задание 5. Тут сравним нашу собственную реализацию случайного леса с sklearn-овской. Для этого воспользуйтесь RandomForestClassifier(n_jobs=1, random_state=17), укажите все те же значения max_depth и max_features, что и раньше. Какое среднее значение ROC AUC на кросс-валидации мы в итоге получили? Выберите самое близкое значение.

Варианты ответа:

  • 0.823
  • 0.833
  • 0.843
  • 0.853

In [22]:
from sklearn.ensemble import RandomForestClassifier

In [23]:
cv_aucs = cross_val_score(RandomForestClassifier(n_estimators=10, max_depth=7, 
                                               max_features=6,
                                               random_state=17, n_jobs=-1,
                                              class_weight='balanced'), 
                        X.values, y.values, scoring="roc_auc", cv=skf)
print("Средняя ROC AUC для случайного леса Sklearn:", np.mean(cv_aucs))


Средняя ROC AUC для случайного леса Sklearn: 0.829143620746

Случайный лес sklearn, подбор гиперпараметров

Задание 6. В 3 задании мы находили оптимальные гиперпараметры для одного дерева, но может быть, для ансамбля эти параметры дерева не будут оптимальными. Давайте проверим это с помощью GridSearchCV (RandomForestClassifier(random_state=17)). Только теперь расширим перебираемые занчения max_depth до 15 включительно, так как в лесу нужны деревья поглубже (а почему именно – вы поняли из статьи). Какими теперь стали лучшие значения гиперпараметров?

Варианты ответа:

  • max_depth=8, max_features=4
  • max_depth=9, max_features=5
  • max_depth=10, max_features=6
  • max_depth=11, max_features=7

In [24]:
max_depth_values = range(5, 15)
max_features_values = [4, 5, 6, 7]
forest_params = {'max_depth': max_depth_values,
               'max_features': max_features_values}

In [25]:
rf = RandomForestClassifier(random_state=17, n_jobs=-1, 
                            class_weight='balanced')
rf_grid_search = GridSearchCV(rf, forest_params, n_jobs=-1, 
                              scoring='roc_auc', cv=skf)
rf_grid_search.fit(X.values, y.values)


Out[25]:
GridSearchCV(cv=StratifiedKFold(n_splits=5, random_state=17, shuffle=True),
       error_score='raise',
       estimator=RandomForestClassifier(bootstrap=True, class_weight='balanced',
            criterion='gini', max_depth=None, max_features='auto',
            max_leaf_nodes=None, min_impurity_decrease=0.0,
            min_impurity_split=None, min_samples_leaf=1,
            min_samples_split=2, min_weight_fraction_leaf=0.0,
            n_estimators=10, n_jobs=-1, oob_score=False, random_state=17,
            verbose=0, warm_start=False),
       fit_params=None, iid=True, n_jobs=-1,
       param_grid={'max_depth': range(5, 15), 'max_features': [4, 5, 6, 7]},
       pre_dispatch='2*n_jobs', refit=True, return_train_score=True,
       scoring='roc_auc', verbose=0)

In [26]:
rf_grid_search.best_score_


Out[26]:
0.83094043147245711

In [27]:
rf_grid_search.best_params_


Out[27]:
{'max_depth': 8, 'max_features': 4}

In [28]:
rf_grid_search.cv_results_["std_test_score"][np.argmax(rf_grid_search.cv_results_["mean_test_score"])]


Out[28]:
0.0023833858843088397

Логистическая регрессия, подбор гиперпараметров

Задание 7. Теперь сравним с логистической регрессией (укажем class_weight='balanced' и random_state=17). Сделайте полный перебор по параметру C из широкого диапазона значений np.logspace(-8, 8, 17). Только сделаем это корректно и выстроим пайплайн – сначала масштабирование, затем обучение модели.

Разберитесь с пайплайнами и проведите кросс-валидацию. Какое получилось лучшее значение средней ROC AUC? Выберите самое близкое значение.

Варианты ответа:

  • 0.778
  • 0.788
  • 0.798
  • 0.808

In [29]:
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression

scaler = StandardScaler()
logit = LogisticRegression(random_state=17, class_weight='balanced')

logit_pipe = Pipeline([('scaler', scaler), ('logit', logit)])
logit_pipe_params = {'logit__C': np.logspace(-8, 8, 17)}

In [30]:
logit_pipe_grid_search = GridSearchCV(logit_pipe, logit_pipe_params, n_jobs=-1, 
                           scoring ='roc_auc', cv=skf)
logit_pipe_grid_search.fit(X.values, y.values)


Out[30]:
GridSearchCV(cv=StratifiedKFold(n_splits=5, random_state=17, shuffle=True),
       error_score='raise',
       estimator=Pipeline(memory=None,
     steps=[('scaler', StandardScaler(copy=True, with_mean=True, with_std=True)), ('logit', LogisticRegression(C=1.0, class_weight='balanced', dual=False,
          fit_intercept=True, intercept_scaling=1, max_iter=100,
          multi_class='ovr', n_jobs=1, penalty='l2', random_state=17,
          solver='liblinear', tol=0.0001, verbose=0, warm_start=False))]),
       fit_params=None, iid=True, n_jobs=-1,
       param_grid={'logit__C': array([  1.00000e-08,   1.00000e-07,   1.00000e-06,   1.00000e-05,
         1.00000e-04,   1.00000e-03,   1.00000e-02,   1.00000e-01,
         1.00000e+00,   1.00000e+01,   1.00000e+02,   1.00000e+03,
         1.00000e+04,   1.00000e+05,   1.00000e+06,   1.00000e+07,
         1.00000e+08])},
       pre_dispatch='2*n_jobs', refit=True, return_train_score=True,
       scoring='roc_auc', verbose=0)

In [31]:
logit_pipe_grid_search.best_score_


Out[31]:
0.78786267573072055

In [32]:
logit_pipe_grid_search.best_params_


Out[32]:
{'logit__C': 100000.0}

In [33]:
logit_pipe_grid_search.cv_results_["std_test_score"][np.argmax(logit_pipe_grid_search.cv_results_["mean_test_score"])]


Out[33]:
0.0048531460319370199

Логистическая регрессия и случайный лес на разреженных признаках

В случае небольшого числа признаков случайный лес показал себя лучше логистической регрессии. Однако один из главных недостатков деревьев проявляется при работе с разреженным данными, например с текстами. Давайте сравним логистическую регрессию и случайный лес в новой задаче. Скачайте данные с отзывами к фильмам отсюда.


In [34]:
# Загрузим данные
df = pd.read_csv("../../data/movie_reviews_train.csv", nrows=50000)

# Разделим данные на текст и целевой признак
X_text = df["text"]
y_text = df["label"]

# Соотношения классов
df.label.value_counts()


Out[34]:
1    32492
0    17508
Name: label, dtype: int64

In [35]:
from sklearn.feature_extraction.text import CountVectorizer
from sklearn.pipeline import Pipeline

# будем разбивать на 3 фолда
skf = StratifiedKFold(n_splits=3, shuffle=True, random_state=17)

# в Pipeline будем сразу преобразовать наш текст и обучать логистическую регрессию
classifier = Pipeline([
    ('vectorizer', CountVectorizer(max_features = 100000, ngram_range = (1, 3))),
    ('clf', LogisticRegression(random_state=17))])

Задание 8. Сделайте полный перебор по параметру C из выборки [0.1, 1, 10, 100]. Какое лучшее значение ROC AUC получилось на кросс-валидации? Выберите самое близкое значение.

Варианты ответа:

  • 0.74
  • 0.75
  • 0.84
  • 0.85

In [36]:
%%time
parameters = {'clf__C': (0.1, 1, 10, 100)}
grid_search = GridSearchCV(classifier, parameters, n_jobs=-1, scoring ='roc_auc', cv=skf)
grid_search = grid_search.fit(X_text, y_text)


CPU times: user 12.6 s, sys: 310 ms, total: 12.9 s
Wall time: 50.6 s

In [37]:
grid_search.best_params_


Out[37]:
{'clf__C': 1}

In [38]:
grid_search.best_score_


Out[38]:
0.85869298711881148

Задание 9. Теперь попробуем сравнить со случайным лесом. Аналогично делаем перебор и получаем максимальное ROC AUC. Выберите самое близкое значение.

Варианты ответа:

  • 0.74
  • 0.75
  • 0.84
  • 0.85

In [39]:
classifier = Pipeline([
    ('vectorizer', CountVectorizer(max_features = 100000, ngram_range = (1, 3))),
    ('clf', RandomForestClassifier(random_state=17, n_jobs=-1))])

min_samples_leaf = [1, 2, 3]
max_features = [0.3, 0.5, 0.7]
max_depth = [None]

In [40]:
%%time
parameters = {'clf__max_features': max_features,
              'clf__min_samples_leaf': min_samples_leaf,
              'clf__max_depth': max_depth}
grid_search = GridSearchCV(classifier, parameters, n_jobs=-1, scoring ='roc_auc', cv=skf)
grid_search = grid_search.fit(X_text, y_text)


CPU times: user 4min 45s, sys: 290 ms, total: 4min 45s
Wall time: 20min 18s

In [41]:
grid_search.best_params_


Out[41]:
{'clf__max_depth': None, 'clf__max_features': 0.5, 'clf__min_samples_leaf': 1}

In [42]:
grid_search.best_score_


Out[42]:
0.74727464051458781